module fin_green_3d_m

    use fin_grn_ext_subs
    implicit none
    private

    public :: FINGREEN3D

    real(rk), parameter :: g = 9.807_rk
    real(rk), parameter :: pi4 = 4*pi
    real(rk), parameter :: ga = 0.5772156649015328_rk
    real(rk), parameter :: eps = 1.0e-6_rk
    complex(rk), parameter :: ci = (0.0_rk, 1.0_rk)

contains

    !>### Description:
    !>
    !> This program computes the three-dimensional free-surface
    !> green function for finite depth and the derivatives of it
    !>
    !> Code original author: yingyi liu       created on  2013.09.07
    !>
    !>### License:
    !>
    !> This routine is part of FINGREEN3D.
    !>
    !> FINGREEN3D is a free software package: you can redistribute it and/or modify it
    !> under the terms of the gnu lesser general public license as published by the
    !> free software foundation, either version 3 of the license, or (at your option)
    !> any later version.
    !>
    !> You should have received a copy of the gnu general public license (either v3
    !> or later), along with FINGREEN3D. if not, see <http://www.gnu.org/licenses/>.
    !>
    !>### Modified on:
    !>
    !> 2022-07-28
    !>
    !>### Reference:
    !>
    !> Yingyi Liu et al. a reliable open-source package for performance evaluation of
    !> floating renewable energy systems in coastal and offshore regions. energy conversion
    !> and management, 2018
    !>
    !> Remark: <br>
    !> The green function can be decomposed into following form,
    !> where the rankine part can be evaluated separately
    !> based on the method of hess and smith (1964) or newman (1986)
    !> in constant panel method due to its high singularity <br>
    !> \( g= 1/r+1/r'+gw \), <br>
    !> where <br>
    !>         \( r  \): \( dsqrt((x-x')^2+(y-y')^2+(z-z')^2) \); <br>
    !>         \( r' \): \( dsqrt((x-x')^2+(y-y')^2+(z+z')^2) \).
    !>
    !>### Parameters:
    !>```text
    !>    input:      rr  --- horizontal distance between the field and the source point
    !>                xp  --- x', first coordinate of the source point
    !>                yf  --- y,  second coordinate of the field point
    !>                yp  --- y', second coordinate of the source point
    !>                zf  --- z,  vertical coordinate of the field point
    !>                zp  --- z', vertical coordinate of the source point
    !>                 v  --- wave number in deep water,i.e., w^2/g
    !>               wvn  --- the first elememnt is the positive root of the dispersion
    !>                        equation v = k*tanh(k*h), the rest elements are the real
    !>                        roots of the equation um*tahn(um*h)= -v
    !>                        (w:circular frequency,g:gravity acceleration)
    !>                nk  --- number of elements in the array wvn
    !>                h   --- water depth (h>0)
    !>                tag --- a flag to exclude the rankine terms or not
    !>                        tag=1, 1/r+1/r1 is excluded;
    !>                        tag=2, 1/r is excluded;
    !>                        tag=3, nothing excluded.
    !>   
    !>    output:    grn(1)  --- value of green's function
    !>               grn(2)  --- derivative of green's function with respect to r
    !>               grn(3)  --- derivative of green's function with respect to z
    !>```
    !>### Contributors list:
    !>
    !> - Yingyi Liu <br>
    !> - 左志华(zuo.zhihua@qq.com)
    pure subroutine FINGREEN3D(rr, zf, zp, v, wvn, nk, h, grn, tag)
        real(rk), intent(in) :: rr          !! Horizontal distance between the field and the source point
        real(rk), intent(in) :: zf          !! Vertical coordinate of the field point
        real(rk), intent(in) :: zp          !! Vertical coordinate of the source point
        real(rk), intent(in) :: v           !! Wave number in deep water,i.e., \( w^2/g \)
        real(rk), intent(in) :: wvn(1:nk)   !! The first elememnt is the positive root of the dispersion equation \( v = k*tanh(k*h) \),
                                            !! The rest elements are the real roots of the equation:
                                            !! \( um*tahn(um*h)= -v \) (w:circular frequency,g:gravity acceleration)
        integer, intent(in) :: nk           !! Number of elements in the array wvn
        real(rk), intent(in) :: h           !! Water depth \( h>0 \)
        integer, intent(in) :: tag          !! A flag to exclude the rankine terms or not <br>
                                            !! tag=1, \( 1/r+1/r1 \) is excluded; <br>
                                            !! tag=2, \( 1/r \) is excluded; <br>
                                            !! tag=3, nothing excluded.
        complex(rk), intent(out) :: grn(1:3)!! Value of green's function <br>
                                            !! grn(1)  --- value of green's function <br>
                                            !! grn(2)  --- derivative of green's function with respect to r <br>
                                            !! grn(3)  --- derivative of green's function with respect to z
        real(rk) r, r1, zpq, rhp, wk
        complex(rk) fg(3)

! initialize parameters based on the inputs
!
        zpq = zf + zp

        r = dsqrt(rr**2 + (zp - zf)**2)
        r1 = dsqrt(rr**2 + zpq*zpq)
        wk = wvn(1)

! choose an appropriate method to calculate the green function
! based on the value of rhp (r/h)
!
        rhp = dabs(rr/h)
        if (rhp < 0.0005_rk) then
            call linton(rr, zf, zp, v, wk, wvn, nk, h, fg)
        elseif (rhp >= 0.0005_rk .and. rhp < 0.05_rk) then
            call pidcock(rr, zf, zp, v, wk, wvn, nk, h, fg)
        elseif (rhp >= 0.05_rk .and. rhp < 0.5_rk) then
            call eigene(rr, zf, zp, wk, wvn, nk, h, fg)
        else
            call eigen(rr, zf, zp, wk, wvn, nk, h, fg)
        end if

! post-process the green function's value
! if r=0, the value of its derivative with respect to r should be zero
!
        grn(:) = dcmplx(0.0_rk, 0.0_rk)

        if (rr /= 0.0_rk) then
            grn(2) = fg(2)
        else
            grn(2) = dcmplx(0.0_rk, 0.0_rk)
        end if

        grn(1) = fg(1)
        grn(3) = fg(3)

! determine whether to exclude the singular terms or not
! this depends on the descritization scheme used in the numerical method
! choose tag=2 or 3 if you use constant panel method
! choose tag=1 if you use higher-order boundary element method
!
        if (tag == 1) then
            continue
        elseif (tag == 2) then
            grn(1) = grn(1) + 1.0_rk/r1
            grn(2) = grn(2) - rr/r1**3
            grn(3) = grn(3) - zpq/r1**3
        elseif (tag == 3) then
            grn(1) = grn(1) + (1.0_rk/r1 + 1.0_rk/r)
            grn(2) = grn(2) - (rr/r1**3 + rr/r**3)
            grn(3) = grn(3) - (zpq/r1**3 + (zf - zp)/r**3)
        end if

    end subroutine FINGREEN3D

!  purpose: this program computes green function and its derivatives
!                using series expansion, without epsilon algorithm, in the region
!                r/h>0.5. the terms needed is 2~10 to 10^-6 accuracy,
!                with a maximum number of term 10.
!
!  license:
!
!    this routine is part of FINGREEN3D.
!
!    FINGREEN3D is a free software package: you can redistribute it and/or modify it
!    under the terms of the gnu lesser general public license as published by the
!    free software foundation, either version 3 of the license, or (at your option)
!    any later version.
!
!    you should have received a copy of the gnu general public license (either v3
!    or later), along with FINGREEN3D. if not, see <http://www.gnu.org/licenses/>.
!
!  code original author:
!
!    yingyi liu on sept 7, 2013
!
!  modified on:
!
!    july 21, 2016
!
!  reference:
!
!    yingyi liu et al. a reliable open-source package for performance evaluation of
!    floating renewable energy systems in coastal and offshore regions. energy conversion
!    and management, 2018
!
!  parameters:
!       input:   r --- horizontal distance of the source point and the field point
!                zf,zp --- z coordinates of the source point and the field point
!                v --- wave number in deep water , i.e., v= w^2/g
!                wk --- positive root of the dispersion equation k*tahn(k*h)= v
!                wvn --- the first elememnt is wk, the rest elements are
!                        the real roots of the equation um*tahn(um*h)= -v
!                nk --- number of elements in the array wvn
!                h --- water depth (h>0)
!       output:  grn --- green function value and its derivatives
!
!  contributors list:
!         yingyi liu
!         to be continued...
!

    pure subroutine eigen(r, zf, zp, wk, wvn, nk, h, grn)
        real(rk), intent(in) :: r, zf, zp, wk, h, wvn(1:nk)
        complex(rk), intent(out) :: grn(3)
        integer, intent(in) :: nk
        integer i, nt
        real(rk) j0, j1, y0, y1, k0, k1
        real(rk) km, pm, rm(3)
        real(rk) sr, sr1, dsrr, dsrz, dsr1r, dsr1z
        complex(rk) p0, dp0, nm0

! initialize some important parameters
!
        nt = 10

        nm0 = h*(1.0_rk + dsinh(2.0_rk*wk*h)/(2.0_rk*wk*h))/2.0_rk
        p0 = pi*dcosh(wk*(zf + h))*dcosh(wk*(zp + h))/nm0
        dp0 = pi*dsinh(wk*(zf + h))*dcosh(wk*(zp + h))/nm0

! calculate the imaginary part of the green function
!
        call jy01b(wk*r, j0, j1, y0, y1)

        grn(1) = p0*(ci*j0 - y0)
        grn(2) = -wk*p0*(ci*j1 - y1)
        grn(3) = wk*dp0*(ci*j0 - y0)

! calculate the real part of the green function
! the trunction terms number depends on the value of rhp
!
        do i = 2, nt

            km = wvn(i)

            pm = 4._rk/((1.0_rk + dsin(2.0_rk*km*h)/(2.0_rk*km*h))*h)

            call ik01a(km*r, k0, k1)

            rm(1) = pm*k0*dcos(km*(zf + h))*dcos(km*(zp + h))
            rm(2) = -pm*km*k1*dcos(km*(zf + h))*dcos(km*(zp + h))
            rm(3) = -pm*km*k0*dsin(km*(zf + h))*dcos(km*(zp + h))

            grn(1) = grn(1) + rm(1)
            grn(2) = grn(2) + rm(2)
            grn(3) = grn(3) + rm(3)

            if (dabs(rm(1)) < eps .and. dabs(rm(2)) < eps .and. dabs(rm(3)) < eps) then
                exit
            end if

        end do

! exclude the singular terms for the purpose of analytical integration over panels
!
        sr = dsqrt(r**2 + (zf - zp)**2)
        sr1 = dsqrt(r**2 + (zf + zp)**2)

        dsrr = r/sr
        dsrz = (zf - zp)/sr
        dsr1r = r/sr1
        dsr1z = (zf + zp)/sr1

        grn(1) = grn(1) - (1.0_rk/sr + 1.0_rk/sr1)
        grn(2) = grn(2) + (dsrr/sr**2 + dsr1r/sr1**2)
        grn(3) = grn(3) + (dsrz/sr**2 + dsr1z/sr1**2)

    end subroutine eigen

!  purpose: this program computes green function and its derivatives
!                using series expansion, with epsilon algorithm, in the region
!                r/h>0.01 and r/h<0.5. the terms needed is 10~80 to
!                10^-5 accuracy, with a maximum number of term 80.
!
!  license:
!
!    this routine is part of FINGREEN3D.
!
!    FINGREEN3D is a free software package: you can redistribute it and/or modify it
!    under the terms of the gnu lesser general public license as published by the
!    free software foundation, either version 3 of the license, or (at your option)
!    any later version.
!
!    you should have received a copy of the gnu general public license (either v3
!    or later), along with FINGREEN3D. if not, see <http://www.gnu.org/licenses/>.
!
!  code original author:
!
!    yingyi liu on sept 7, 2013
!
!  modified on:
!
!    july 21, 2016
!
!  reference:
!
!    yingyi liu et al. a reliable open-source package for performance evaluation of
!    floating renewable energy systems in coastal and offshore regions. energy conversion
!    and management, 2018
!
!  parameters:
!       input:   r --- horizontal distance of the source point and the field point
!                zf,zp --- z coordinates of the source point and the field point
!                v --- wave number in deep water , i.e., v= w^2/g
!                wk --- positive root of the dispersion equation k*tahn(k*h)= v
!                wvn --- the first elememnt is wk, the rest elements are
!                        the real roots of the equation um*tahn(um*h)= -v
!                nk --- number of elements in the array wvn
!                h --- water depth (h>0)
!       output:  grn --- green function value and its derivatives
!
!  contributors list:
!         yingyi liu
!         to be continued...
!

    pure subroutine eigene(r, zf, zp, wk, wvn, nk, h, grn)
        real(rk), intent(in) :: r, zf, zp, h, wk, wvn(1:nk)
        integer, intent(in) :: nk
        complex(rk), intent(out) :: grn(3)

        integer i, nt, i_pade, k_pade
        real(rk) j0, j1, y0, y1, k0, k1
        real(rk) km, pm, rm(3), rl(3), err, rhp
        real(rk) sr, sr1, dsrr, dsrz, dsr1r, dsr1z
        real(rk) g(0:nk - 1), gr(0:nk - 1), gz(0:nk - 1)
        complex(rk) p0, dp0, nm0

! initialize some important parameters
!
        rhp = dabs(r/h)
        nt = int(-88.89*rhp + 54.45)

        nm0 = h*(1.0_rk + dsinh(2.0_rk*wk*h)/(2.0_rk*wk*h))/2.0_rk
        p0 = pi*dcosh(wk*(zf + h))*dcosh(wk*(zp + h))/nm0
        dp0 = pi*dsinh(wk*(zf + h))*dcosh(wk*(zp + h))/nm0

! calculate the imaginary part of the green function
!
        call jy01b(wk*r, j0, j1, y0, y1)

        grn(1) = p0*(ci*j0 - y0)
        grn(2) = -wk*p0*(ci*j1 - y1)
        grn(3) = wk*dp0*(ci*j0 - y0)

! calculate the real part of the green function
! the trunction terms number depends on the value of rhp
!
        g(0) = 0.0_rk
        gr(0) = 0.0_rk
        gz(0) = 0.0_rk

        do i = 2, nt

            km = wvn(i)

            pm = 4._rk/((1.0_rk + dsin(2.0_rk*km*h)/(2.0_rk*km*h))*h)

            call ik01a(km*r, k0, k1)

            rm(1) = pm*k0*dcos(km*(zf + h))*dcos(km*(zp + h))
            rm(2) = -pm*km*k1*dcos(km*(zf + h))*dcos(km*(zp + h))
            rm(3) = -pm*km*k0*dsin(km*(zf + h))*dcos(km*(zp + h))

            g(i - 1) = g(i - 2) + rm(1)
            gr(i - 1) = gr(i - 2) + rm(2)
            gz(i - 1) = gz(i - 2) + rm(3)

        end do

! call the epsilon algorithm to approximate the limit of the infinite series
!
        call limes(nt - 1, g(0:nt - 1), rl(1), i_pade, k_pade, err)
        call limes(nt - 1, gr(0:nt - 1), rl(2), i_pade, k_pade, err)
        call limes(nt - 1, gz(0:nt - 1), rl(3), i_pade, k_pade, err)

        grn(1) = grn(1) + rl(1)
        grn(2) = grn(2) + rl(2)
        grn(3) = grn(3) + rl(3)

! exclude the singular terms for the purpose of analytical integration over panels
!
        sr = dsqrt(r**2 + (zf - zp)**2)
        sr1 = dsqrt(r**2 + (zf + zp)**2)

        dsrr = r/sr
        dsrz = (zf - zp)/sr
        dsr1r = r/sr1
        dsr1z = (zf + zp)/sr1

        grn(1) = grn(1) - (1.0_rk/sr + 1.0_rk/sr1)
        grn(2) = grn(2) + (dsrr/sr**2 + dsr1r/sr1**2)
        grn(3) = grn(3) + (dsrz/sr**2 + dsr1z/sr1**2)

    end subroutine eigene

!  purpose: this program computes green function and its derivatives
!                using pidcock's method, with epsilon algorithm, in the region
!                r/h>0.0005 and r/h<0.01. the terms needed is 50~100 to
!                10^-4 accuracy, with a maximum number of term 100.
!
!  license:
!
!    this routine is part of FINGREEN3D.
!
!    FINGREEN3D is a free software package: you can redistribute it and/or modify it
!    under the terms of the gnu lesser general public license as published by the
!    free software foundation, either version 3 of the license, or (at your option)
!    any later version.
!
!    you should have received a copy of the gnu general public license (either v3
!    or later), along with FINGREEN3D. if not, see <http://www.gnu.org/licenses/>.
!
!  code original author:
!
!    yingyi liu on sept 7, 2013
!
!  modified on:
!
!    july 21, 2016
!
!  reference:
!
!    yingyi liu et al. a reliable open-source package for performance evaluation of
!    floating renewable energy systems in coastal and offshore regions. energy conversion
!    and management, 2018
!
!  parameters:
!       input:   r --- horizontal distance of the source point and the field point
!                zf,zp --- z coordinates of the source point and the field point
!                v --- wave number in deep water , i.e., v= w^2/g
!                wk --- positive root of the dispersion equation k*tahn(k*h)= v
!                wvn --- the first elememnt is wk, the rest elements are
!                        the real roots of the equation um*tahn(um*h)= -v
!                nk --- number of elements in the array wvn
!                h --- water depth (h>0)
!       output:  grn --- green function value and its derivatives
!
!  contributors list:
!         yingyi liu
!         to be continued...
!

    pure subroutine pidcock(r, zf, zp, v, wk, wvn, nk, h, grn)
        real(rk), intent(in) :: r, zf, zp, h, v, wk, wvn(1:nk)
        integer, intent(in) :: nk
        complex(rk), intent(out) :: grn(3)

        integer nt, i, i_pade, k_pade

        real(rk) j0, j1, y0, y1, k0, k0t, k1, k1t
        real(rk) nm0, p0, dp0, x
        real(rk) rm(3), sm(3), rl(3), err, rhp
        real(rk) um, umt, mp, ac, am, amt, dam, damt
        real(rk) sr, sr1, dsrr, dsrz, dsr1r, dsr1z
        real(rk) g(0:nk - 1), gr(0:nk - 1), gz(0:nk - 1)

! initialize some important parameters
!
        nm0 = h*(1.0_rk + dsinh(2.0_rk*wk*h)/(2.0_rk*wk*h))/2.0_rk
        p0 = pi*dcosh(wk*(zf + h))*dcosh(wk*(zp + h))/nm0
        dp0 = wk*pi*dsinh(wk*(zf + h))*dcosh(wk*(zp + h))/nm0

! calculate the imaginary part of the green function
!
        call jy01b(wk*r, j0, j1, y0, y1)

        grn(1) = p0*(ci*j0 - y0) + 2.0_rk/h*(ga + dlog(r/4._rk/h))
        grn(2) = -wk*p0*(ci*j1 - y1) + 2.0_rk/(r*h)
        grn(3) = dp0*(ci*j0 - y0)

! calculate the real part, for the infinite series of eq.(13)
! the series applies an acceleration method by substracting an asymptotic series
!
        rhp = dabs(r/h)
        nt = int(-1010.10*rhp + 100.50)

        g(0) = 0.0_rk
        gr(0) = 0.0_rk
        gz(0) = 0.0_rk

        do i = 2, nt

            mp = (i - 1)*pi
            umt = mp/h
            um = wvn(i)
            x = um*r

            call ik01a(x, k0, k1)
            call ik01a(umt*r, k0t, k1t)

            ac = (um**2 + v**2)/((um**2 + v**2)*h - v)
            am = ac*dcos(um*(zf + h))*dcos(um*(zp + h))
            amt = dcos(umt*(zf + h))*dcos(umt*(zp + h))/h
            dam = -um*ac*dsin(um*(zf + h))*dcos(um*(zp + h))
            damt = -umt*dsin(umt*(zf + h))*dcos(umt*(zp + h))/h

            rm(1) = am*k0 - amt*k0t
            rm(2) = -um*am*k1 + umt*amt*k1t
            rm(3) = dam*k0 - damt*k0t

            g(i - 1) = g(i - 2) + rm(1)
            gr(i - 1) = gr(i - 2) + rm(2)
            gz(i - 1) = gz(i - 2) + rm(3)

        end do

! call the epsilon algorithm to approximate the limit of the infinite series
!
        call limes(nt - 1, g(0:nt - 1), rl(1), i_pade, k_pade, err)
        call limes(nt - 1, gr(0:nt - 1), rl(2), i_pade, k_pade, err)
        call limes(nt - 1, gz(0:nt - 1), rl(3), i_pade, k_pade, err)

        sm(1) = 4._rk*rl(1)
        sm(2) = 4._rk*rl(2)
        sm(3) = 4._rk*rl(3)

! calculate the sum of rankine terms by chebyshev approximations
!
        grn(1) = grn(1) + sm(1) + imgs(r, zf, zp, h)
        grn(2) = grn(2) + sm(2) + dgsr(r, zf, zp, h)
        grn(3) = grn(3) + sm(3) + dgsz(r, zf, zp, h)

! exclude the singular terms for the purpose of analytical integration over panels
!
        sr = dsqrt(r**2 + (zf - zp)**2)
        sr1 = dsqrt(r**2 + (zf + zp)**2)

        dsrr = r/sr
        dsrz = (zf - zp)/sr
        dsr1r = r/sr1
        dsr1z = (zf + zp)/sr1

        grn(1) = grn(1) - (1.0_rk/sr + 1.0_rk/sr1)
        grn(2) = grn(2) + (dsrr/sr**2 + dsr1r/sr1**2)
        grn(3) = grn(3) + (dsrz/sr**2 + dsr1z/sr1**2)

    end subroutine pidcock

!  purpose: this program computes green function and its derivatives
!                using linton's method, with epsilon algorithm, in the region
!                r/h<0.0005. the terms needed is about 3~10 to
!                10^-6 accuracy, with a maximum number of term 10.
!
!  license:
!
!    this routine is part of FINGREEN3D.
!
!    FINGREEN3D is a free software package: you can redistribute it and/or modify it
!    under the terms of the gnu lesser general public license as published by the
!    free software foundation, either version 3 of the license, or (at your option)
!    any later version.
!
!    you should have received a copy of the gnu general public license (either v3
!    or later), along with FINGREEN3D. if not, see <http://www.gnu.org/licenses/>.
!
!  code original author:
!
!    yingyi liu on sept 7, 2013
!
!  modified on:
!
!    july 21, 2016
!
!  reference:
!
!    yingyi liu et al. a reliable open-source package for performance evaluation of
!    floating renewable energy systems in coastal and offshore regions. energy conversion
!    and management, 2018
!
!  parameters:
!       input:   r --- horizontal distance of the source point and the field point
!                zf,zp --- z coordinates of the source point and the field point
!                v --- wave number in deep water , i.e., v= w^2/g
!                wk --- positive root of the dispersion equation k*tahn(k*h)= v
!                wvno --- the first elememnt is wk, the rest elements are
!                        the real roots of the equation um*tahn(um*h)= -v
!                nk --- number of elements in the array wvno
!                h --- water depth (h>0)
!       output:  grn --- green function value and its derivatives
!
!  contributors list:
!         yingyi liu
!         to be continued...
!

    pure subroutine linton(r, zf, zp, v, wk, wvno, nk, h, grn)
        real(rk), intent(in) :: r, zf, zp, h, v, wk, wvno(nk)
        integer, intent(in) :: nk
        complex(rk), intent(out) :: grn(3)

        integer nt, i, flag
        real(rk) pa, a
        real(rk) sr, sr1, sr2, dsrr, dsrz, dsr1r, dsr1z, dsr2r, dsr2z
        real(rk) ro(4), dror(4), droz(4), errsum(3), bj0, bj1
        real(rk) cof(nk), dcof(nk)
        complex(rk) km, rm(3)

!-----------------------------------------------------------------------------
!  i. optimal selection of the parameter a
!-----------------------------------------------------------------------------

        flag = 2

        if (flag == 1) then
            a = 0.25_rk
        elseif (flag == 2) then
            if (v < 0.02_rk) then
                a = 3324.2_rk*v**3 - 258.89_rk*v**2 + 8.3958_rk*v + 1.3208_rk
            elseif (v < 0.1_rk) then
                a = -30.303_rk*v**3 - 2.8138_rk*v**2 + 2.7740_rk*v + 0.05052_rk
            elseif (v < 1.0_rk) then
                a = 1.5699_rk*v**3 - 2.5075_rk*v**2 + 1.8316_rk*v + 0.1077_rk
            else
                a = -0.000069641_rk*v**3 - 0.0017082_rk*v**2 + 0.45566_rk*v + 0.80551_rk
            end if
            a = dabs(a)/(wk*h)
        end if
!
!       print*,a
!
!-----------------------------------------------------------------------------
!  ii. initialisation of various parameters
!-----------------------------------------------------------------------------

        pa = a*h
        sr = dsqrt(r**2 + (zf - zp)**2)
        sr1 = dsqrt(r**2 + (2.0_rk*h + zf + zp)**2)
        sr2 = dsqrt(r**2 + (zf + zp)**2)

        call rchi(r, zf, zp, h, ro, dror, droz)

        dsrr = r/sr
        dsrz = (zf - zp)/sr
        dsr1r = r/sr1
        dsr1z = (2.0_rk*h + zf + zp)/sr1
        dsr2r = r/sr2
        dsr2z = (zf + zp)/sr2
!
!-----------------------------------------------------------------------------
!  iii. calculating the part of evanescent modes
!-----------------------------------------------------------------------------

        grn(:) = (0.0_rk, 0.0_rk)
        nt = 20

        call coef(r, pa, a, nt, wvno, wk, cof)
        call dcoef(r, pa, a, nt, wvno, wk, dcof)

        do i = 1, nt

            if (i == 1) then
                km = -ci*wk
            else
                km = wvno(i)
            end if

            rm(1) = -cof(i)*cdcos(km*(zf + h))*cdcos(km*(zp + h))/nm(h, km)
            rm(2) = -dcof(i)*cdcos(km*(zf + h))*cdcos(km*(zp + h))/nm(h, km)
            rm(3) = cof(i)*km*cdsin(km*(zf + h))*cdcos(km*(zp + h))/nm(h, km)

            grn(1) = grn(1) + rm(1)
            grn(2) = grn(2) + rm(2)
            grn(3) = grn(3) + rm(3)

            if (i > 2 .and. cdabs(rm(1)) < eps .and. cdabs(rm(2)) < eps .and. cdabs(rm(3)) < eps) then
                exit
            end if

        end do
!
!-----------------------------------------------------------------------------
!  iv. calculating the part of error functions
!-----------------------------------------------------------------------------

        errsum(:) = 0.0_rk

        do i = 1, 4

            errsum(1) = errsum(1) + erfcc(ro(i)/pa)/(pi4*ro(i))
            errsum(2) = errsum(2) + (derfcc(ro(i)/pa)*dror(i)*ro(i)/pa - dror(i)*erfcc(ro(i)/pa))/(pi4*ro(i)**2)
            errsum(3) = errsum(3) + (derfcc(ro(i)/pa)*droz(i)*ro(i)/pa - droz(i)*erfcc(ro(i)/pa))/(pi4*ro(i)**2)

        end do

        if (dabs(sr) > eps) then

            errsum(1) = errsum(1) + (erfcc(sr/pa) - 1.0_rk)/(pi4*sr) + erfcc(sr1/pa)/(pi4*sr1) - 1.0_rk/(pi4*sr2)
            errsum(2) = errsum(2) - (erfcc(sr/pa) - 1.0_rk)/(pi4*sr**2)*dsrr + derfcc(sr/pa)/(pi4*sr*pa)*dsrr - erfcc(sr1/pa) &
                        *dsr1r/(pi4*sr1**2) + derfcc(sr1/pa)*dsr1r/(pi4*sr1*pa) + dsr2r/(pi4*sr2**2)
            errsum(3) = errsum(3) - (erfcc(sr/pa) - 1.0_rk)/(pi4*sr**2)*dsrz + derfcc(sr/pa)/(pi4*sr*pa)*dsrz - erfcc(sr1/pa) &
                        *dsr1z/(pi4*sr1**2) + derfcc(sr1/pa)*dsr1z/(pi4*sr1*pa) + dsr2z/(pi4*sr2**2)

        else

            errsum(1) = errsum(1) + erfcc(sr1/pa)/(pi4*sr1) - 1.0_rk/(pi4*sr2)
            errsum(2) = errsum(2) - erfcc(sr1/pa)*dsr1r/(pi4*sr1**2) + derfcc(sr1/pa)*dsr1r/(pi4*sr1*pa) + dsr2r/(pi4*sr2**2)
            errsum(3) = errsum(3) - erfcc(sr1/pa)*dsr1z/(pi4*sr1**2) + derfcc(sr1/pa)*dsr1z/(pi4*sr1*pa) + dsr2z/(pi4*sr2**2)

        end if
!
!-----------------------------------------------------------------------------
!  v. calculating the part of integration
!-----------------------------------------------------------------------------

        grn(1) = grn(1) - errsum(1) - v*aquadf(v, r, zf, zp, h, 0.0_rk, pa/2.0_rk, eps, 1)/(2.0_rk*pi)
        grn(2) = grn(2) - errsum(2) + v*aquadf(v, r, zf, zp, h, 0.0_rk, pa/2.0_rk, eps, 2)/(2.0_rk*pi)
        grn(3) = grn(3) - errsum(3) - v*aquadf(v, r, zf, zp, h, 0.0_rk, pa/2.0_rk, eps, 3)/(2.0_rk*pi)

        call jy01bj(wk*r, bj0, bj1)

        km = -ci*wk
        grn(1) = grn(1) - ci*bj0*dcosh(wk*(zf + h))*dcosh(wk*(zp + h))/(4._rk*nm(h, km))
        grn(2) = grn(2) + ci*wk*bj1*dcosh(wk*(zf + h))*dcosh(wk*(zp + h))/(4._rk*nm(h, km))
        grn(3) = grn(3) - ci*wk*bj0*dsinh(wk*(zf + h))*dcosh(wk*(zp + h))/(4._rk*nm(h, km))

        grn(1) = -pi4*grn(1)
        grn(2) = -pi4*grn(2)
        grn(3) = -pi4*grn(3)

    end subroutine linton

!  purpose: this program computes rankine image source
!
!  license:
!
!    this routine is part of FINGREEN3D.
!
!    FINGREEN3D is a free software package: you can redistribute it and/or modify it
!    under the terms of the gnu lesser general public license as published by the
!    free software foundation, either version 3 of the license, or (at your option)
!    any later version.
!
!    you should have received a copy of the gnu general public license (either v3
!    or later), along with FINGREEN3D. if not, see <http://www.gnu.org/licenses/>.
!
!  code original author:
!
!    yingyi liu on sept 7, 2013
!
!  modified on:
!
!    july 21, 2016
!
!  reference:
!
!    yingyi liu et al. a reliable open-source package for performance evaluation of
!    floating renewable energy systems in coastal and offshore regions. energy conversion
!    and management, 2018
!
!  parameters:
!       input:   r --- horizontal distance of the source point and the field point
!                zf,zp --- z coordinates of the source point and the field point
!                h --- water depth (h>0)
!       output:  imgs --- value of the rankine image source in pidcock's method
!
!  contributors list:
!         yingyi liu
!         to be continued...
!

    pure real(rk) function imgs(r, zf, zp, h)
        real(rk), intent(in) :: r, zf, zp, h
        integer l, m, n
        real(rk) t1, t2, t3, sum1, sum2

        real(rk), parameter :: amn(0:4, 0:4) = &
                               reshape([0.0_rk, -0.02525711_rk, 0.00086546_rk, -0.00004063_rk, &
                                        0.00000193_rk, 0.05051418_rk, -0.00692380_rk, 0.00073292_rk, &
                                        -0.00006636_rk, 0.00000398_rk, 0.00230838_rk, -0.00097875_rk, &
                                        0.00020597_rk, -0.00003333_rk, 0.00000524_rk, 0.00012934_rk, &
                                        -0.00010879_rk, 0.00003965_rk, -0.00000891_rk, 0.0_rk, &
                                        0.00000913_rk, -0.00001270_rk, 0.00000466_rk, 0.0_rk, 0.0_rk], shape(amn))

        t1 = r/h
        t2 = (zf - zp)/h
        t3 = (zf + zp + 2.0_rk*h)/h

        sum1 = 0.0_rk
        do l = -1, 1
            sum1 = sum1 + 1.0_rk/rs(t1, t2 + 2.0_rk*l) + 1.0_rk/rs(t1, t3 + 2.0_rk*l)
        end do

        sum2 = 0.0_rk
        do m = 0, 4
            do n = 0, 4
                sum2 = sum2 + amn(m, n)*t1**(2.0_rk*m)*(t2**(2.0_rk*n) + t3**(2.0_rk*n))
            end do
        end do

        imgs = (sum1 - 2.0_rk + sum2)/h

    end function imgs

!  purpose: this program computes the derivative of rankine image source
!                with respect to r
!
!  license:
!
!    this routine is part of FINGREEN3D.
!
!    FINGREEN3D is a free software package: you can redistribute it and/or modify it
!    under the terms of the gnu lesser general public license as published by the
!    free software foundation, either version 3 of the license, or (at your option)
!    any later version.
!
!    you should have received a copy of the gnu general public license (either v3
!    or later), along with FINGREEN3D. if not, see <http://www.gnu.org/licenses/>.
!
!  code original author:
!
!    yingyi liu on sept 7, 2013
!
!  modified on:
!
!    july 21, 2016
!
!  reference:
!
!    yingyi liu et al. a reliable open-source package for performance evaluation of
!    floating renewable energy systems in coastal and offshore regions. energy conversion
!    and management, 2018
!
!  parameters:
!       input:   r --- horizontal distance of the source point and the field point
!                zf,zp --- z coordinates of the source point and the field point
!                h --- water depth (h>0)
!       output:  dgsr --- derivative of the rankine image source with respect to r
!
!  contributors list:
!         yingyi liu
!         to be continued...
!

    pure real(rk) function dgsr(r, zf, zp, h)
        real(rk), intent(in) :: r, zf, zp, h
        integer l, m, n
        real(rk) t1, t2, t3, sum1, sum2

        real(rk), parameter :: amn(0:4, 0:4) = &
                               reshape([0.0_rk, -0.02525711_rk, 0.00086546_rk, -0.00004063_rk, &
                                        0.00000193_rk, 0.05051418_rk, -0.00692380_rk, 0.00073292_rk, &
                                        -0.00006636_rk, 0.00000398_rk, 0.00230838_rk, -0.00097875_rk, &
                                        0.00020597_rk, -0.00003333_rk, 0.00000524_rk, 0.00012934_rk, &
                                        -0.00010879_rk, 0.00003965_rk, -0.00000891_rk, 0.0_rk, &
                                        0.00000913_rk, -0.00001270_rk, 0.00000466_rk, 0.0_rk, 0.0_rk], shape(amn))

        t1 = r/h
        t2 = (zf - zp)/h
        t3 = (zf + zp + 2.0_rk*h)/h

        sum1 = 0.0_rk
        do l = -1, 1
            sum1 = sum1 + t1/(rs(t1, t2 + 2.0_rk*l))**3 + t1/(rs(t1, t3 + 2.0_rk*l))**3
        end do

        sum2 = 0.0_rk
        do m = 1, 4
            do n = 0, 4
                sum2 = sum2 + 2.0_rk*m*amn(m, n)*t1**(2.0_rk*m - 1.0_rk)*(t2**(2.0_rk*n) + t3**(2.0_rk*n))
            end do
        end do

        dgsr = (-sum1 + sum2)/h**2

    end function dgsr

!  purpose: this program computes the derivative of rankine image source
!                with respect to z
!
!  license:
!
!    this routine is part of FINGREEN3D.
!
!    FINGREEN3D is a free software package: you can redistribute it and/or modify it
!    under the terms of the gnu lesser general public license as published by the
!    free software foundation, either version 3 of the license, or (at your option)
!    any later version.
!
!    you should have received a copy of the gnu general public license (either v3
!    or later), along with FINGREEN3D. if not, see <http://www.gnu.org/licenses/>.
!
!  code original author:
!
!    yingyi liu on sept 7, 2013
!
!  modified on:
!
!    july 21, 2016
!
!  reference:
!
!    yingyi liu et al. a reliable open-source package for performance evaluation of
!    floating renewable energy systems in coastal and offshore regions. energy conversion
!    and management, 2018
!
!  parameters:
!       input:   r --- horizontal distance of the source point and the field point
!                zf,zp --- z coordinates of the source point and the field point
!                h --- water depth (h>0)
!       output:  dgsz --- derivative of the rankine image source with respect to z
!
!  contributors list:
!         yingyi liu
!         to be continued...
!

    pure real(rk) function dgsz(r, zf, zp, h)
        real(rk), intent(in) :: r, zf, zp, h
        integer l, m, n
        real(rk) t1, t2, t3, sum1, sum2

        real(rk), parameter :: amn(0:4, 0:4) = &
                               reshape([0.0_rk, -0.02525711_rk, 0.00086546_rk, -0.00004063_rk, &
                                        0.00000193_rk, 0.05051418_rk, -0.00692380_rk, 0.00073292_rk, &
                                        -0.00006636_rk, 0.00000398_rk, 0.00230838_rk, -0.00097875_rk, &
                                        0.00020597_rk, -0.00003333_rk, 0.00000524_rk, 0.00012934_rk, &
                                        -0.00010879_rk, 0.00003965_rk, -0.00000891_rk, 0.0_rk, &
                                        0.00000913_rk, -0.00001270_rk, 0.00000466_rk, 0.0_rk, 0.0_rk], shape(amn))

        t1 = r/h
        t2 = (zf - zp)/h
        t3 = (zf + zp + 2.0_rk*h)/h

        sum1 = 0.0_rk
        do l = -1, 1
            sum1 = sum1 + (t2 + 2.0_rk*l)/(rs(t1, t2 + 2.0_rk*l))**3 + (t3 + 2.0_rk*l)/(rs(t1, t3 + 2.0_rk*l))**3
        end do

        sum2 = 0.0_rk
        do m = 0, 4
            do n = 1, 4
                sum2 = sum2 + 2.0_rk*n*amn(m, n)*t1**(2.0_rk*m)*(t2**(2.0_rk*n - 1.0_rk) + t3**(2.0_rk*n - 1.0_rk))
            end do
        end do

        dgsz = (-sum1 + sum2)/h**2

    end function dgsz

!
!    ======================================================
!          function dsqrt(r^2+z^2)
!    ======================================================

    pure real(rk) function rs(r, z)

        real(rk), intent(in) :: r, z

        rs = dsqrt(r**2 + z**2)

    end function rs

!  purpose: this program computes the expansion coefficients in
!                linton's method
!
!  license:
!
!    this routine is part of FINGREEN3D.
!
!    FINGREEN3D is a free software package: you can redistribute it and/or modify it
!    under the terms of the gnu lesser general public license as published by the
!    free software foundation, either version 3 of the license, or (at your option)
!    any later version.
!
!    you should have received a copy of the gnu general public license (either v3
!    or later), along with FINGREEN3D. if not, see <http://www.gnu.org/licenses/>.
!
!  code original author:
!
!    yingyi liu on sept 7, 2013
!
!  modified on:
!
!    july 21, 2016
!
!  reference:
!
!    yingyi liu et al. a reliable open-source package for performance evaluation of
!    floating renewable energy systems in coastal and offshore regions. energy conversion
!    and management, 2018
!
!  parameters:
!       input:   r --- horizontal distance of the source point and the field point
!                pa --- a key parameter determing convergence, i.e., pa= a*h
!                a ---  a key parameter related with convergence
!                m --- number of elements in the array wvno
!                wvno --- the first elememnt is wk, the rest elements are
!                        the real roots of the equation um*tahn(um*h)= -v
!                wk --- positive root of the dispersion equation k*tahn(k*h)= v
!       output:  cof --- expansion coefficients in linton's expression
!
!  contributors list:
!         yingyi liu
!         to be continued...
!

    pure subroutine coef(r, pa, a, m, wvno, wk, cof)
        integer i
        integer, intent(in) :: m
        real(rk), intent(in) :: r, pa, a, wvno(m), wk
        real(rk), intent(out) :: cof(m)
        real(rk) by0, dnt, step, p1, e1, ei
        step = 100._rk
        dnt = pa**2/4._rk

! calculate the first expansion coefficient a0
!
        if (dabs(r) < eps) then

            call eix(wk**2*dnt, ei)
            cof(1) = -ei/pi4

        elseif (r/a <= 0.5_rk) then

            cof(1) = sa0(r, pa, wk)

        else

            call jy01by0(wk*r, by0)
            cof(1) = -by0/4._rk - aquad(r, wk, 1, 0.0_rk, dnt, eps, 1)

        end if

! calculate the rest expansion coefficients ai (i=2,m)
! the truncation terms number depends on whether the series
! converged to a specified accuracy or not

        do i = 2, m

            if (dabs(r) < eps) then

                call e1xa(wvno(i)**2*dnt, e1)
                cof(i) = e1/pi4

            elseif (r/a <= 1.0_rk) then

                cof(i) = sam(r, pa, wvno(i))

            else

                p1 = dnt + step
                do while (dabs(f1(r, wvno(i), i, p1)) > eps)
                    p1 = p1 + step
                end do

                cof(i) = aquad(r, wvno(i), i, dnt, p1, eps, 1)

            end if

        end do

    end subroutine coef

!  purpose: this program computes derivative of the expansion coefficients
!                with repect to r in linton's method
!
!  license:
!
!    this routine is part of FINGREEN3D.
!
!    FINGREEN3D is a free software package: you can redistribute it and/or modify it
!    under the terms of the gnu lesser general public license as published by the
!    free software foundation, either version 3 of the license, or (at your option)
!    any later version.
!
!    you should have received a copy of the gnu general public license (either v3
!    or later), along with FINGREEN3D. if not, see <http://www.gnu.org/licenses/>.
!
!  code original author:
!
!    yingyi liu on sept 7, 2013
!
!  modified on:
!
!    july 21, 2016
!
!  reference:
!
!    yingyi liu et al. a reliable open-source package for performance evaluation of
!    floating renewable energy systems in coastal and offshore regions. energy conversion
!    and management, 2018
!
!  parameters:
!       input:   r --- horizontal distance of the source point and the field point
!                pa --- a key parameter determing convergence, i.e., pa= a*h
!                a ---  a key parameter related with convergence
!                m --- number of elements in the array wvno
!                wvno --- the first elememnt is wk, the rest elements are
!                        the real roots of the equation um*tahn(um*h)= -v
!                wk --- positive root of the dispersion equation k*tahn(k*h)= v
!       output:  dcoef --- derivative of expansion coefficients with repect to r
!
!  contributors list:
!         yingyi liu
!         to be continued...
!

    pure subroutine dcoef(r, pa, a, m, wvno, wk, dcof)
        integer i
        integer, intent(in) :: m
        real(rk), intent(in) :: r, pa, a, wvno(m), wk
        real(rk), intent(out) :: dcof(m)
        real(rk) by1, dnt, step, p1

        step = 100._rk
        dnt = pa**2/4._rk

! calculate the derivative of the first expansion coefficient a0
!
        if (dabs(r) < eps) then

            dcof(1) = 0.0_rk

        elseif (r/a <= 0.5_rk) then

            dcof(1) = dsa0(r, pa, wk)

        else

            call jy01by1(wk*r, by1)
            dcof(1) = wk*by1/4._rk - aquad(r, wk, 1, 0.0_rk, dnt, eps, 2)

        end if

! calculate the derivative of the rest expansion coefficients ai (i=2,m)
! the truncation terms number depends on whether the series
! converged to a specified accuracy or not

        do i = 2, m

            if (dabs(r) < eps) then

                dcof(i) = 0.0_rk

            elseif (r/a <= 1.0_rk) then

                dcof(i) = dsam(r, pa, wvno(i))

            else

                p1 = dnt + step
                do while (dabs(f2(r, wvno(i), i, p1)) > eps)
                    p1 = p1 + step
                end do

                dcof(i) = aquad(r, wvno(i), i, dnt, p1, eps, 2)

            end if

        end do

    end subroutine dcoef

!  purpose: this program computes the series representations of expansion
!                coefficients (m=0)
!
!  license:
!
!    this routine is part of FINGREEN3D.
!
!    FINGREEN3D is a free software package: you can redistribute it and/or modify it
!    under the terms of the gnu lesser general public license as published by the
!    free software foundation, either version 3 of the license, or (at your option)
!    any later version.
!
!    you should have received a copy of the gnu general public license (either v3
!    or later), along with FINGREEN3D. if not, see <http://www.gnu.org/licenses/>.
!
!  code original author:
!
!    yingyi liu on sept 7, 2013
!
!  modified on:
!
!    july 17, 2018
!
!  reference:
!
!    yingyi liu et al. a reliable open-source package for performance evaluation of
!    floating renewable energy systems in coastal and offshore regions. energy conversion
!    and management, 2018
!
!  parameters:
!       input:   r --- horizontal distance of the source point and the field point
!                pa --- a key parameter determing convergence, i.e., pa= a*h
!                a ---  a key parameter related with convergence
!                m --- number of elements in the array wvno
!                wvno --- the first elememnt is wk, the rest elements are
!                        the real roots of the equation um*tahn(um*h)= -v
!                wk --- positive root of the dispersion equation k*tahn(k*h)= v
!       output:  sa0 --- value of first expansion coefficient by series representation (m=0)
!
!  contributors list:
!         yingyi liu
!         to be continued...
!

    pure real(rk) function sa0(r, pa, wk)
        real(rk), intent(in) :: r, pa, wk
        integer n, p
        real(rk) psi, rm, en, sgn, facn

        n = 1
        rm = 100._rk
        sa0 = -ga - 2.0_rk*dlog(wk*pa/2.0_rk)
        sgn = -1.0_rk
        facn = 1.0_rk

! loop: the summation stops at the condition when rm meets the tolerance
!
        do while (dabs(rm) > eps)

            call enxa(n + 1, (r/pa)**2, en)

            psi = 0.0_rk
            do p = 1, n
                psi = psi + 1.0_rk/p
            end do

! calculate the value of each term in the infinite series expansion
!
            rm = sgn*(r/pa)**(2.0_rk*n)/(facn*n) - (wk*pa/2.0_rk)**(2.0_rk*n)*en/facn - 2.0_rk*sgn &
                 *(wk*r/2.0_rk)**(2.0_rk*n)*(dlog(wk*r/2.0_rk) + ga - psi)/facn**2

            !rm=-(wk*pa/2.0_rk)**(2.0_rk*n)/facn*en

            sa0 = sa0 + rm

            n = n + 1
            sgn = -sgn   ! sign function depends on the power n
            facn = facn*n   ! factorial of the numbers from 1 to n

        end do

        sa0 = sa0/pi4

    end function sa0

!  purpose: this program computes the series representations of expansion
!                coefficients with repect to r in linton's method (m=0)
!
!  license:
!
!    this routine is part of FINGREEN3D.
!
!    FINGREEN3D is a free software package: you can redistribute it and/or modify it
!    under the terms of the gnu lesser general public license as published by the
!    free software foundation, either version 3 of the license, or (at your option)
!    any later version.
!
!    you should have received a copy of the gnu general public license (either v3
!    or later), along with FINGREEN3D. if not, see <http://www.gnu.org/licenses/>.
!
!  code original author:
!
!    yingyi liu on sept 7, 2013
!
!  modified on:
!
!    july 17, 2018
!
!  reference:
!
!    yingyi liu et al. a reliable open-source package for performance evaluation of
!    floating renewable energy systems in coastal and offshore regions. energy conversion
!    and management, 2018
!
!  parameters:
!       input:   r --- horizontal distance of the source point and the field point
!                pa --- a key parameter determing convergence, i.e., pa= a*h
!                a ---  a key parameter related with convergence
!                m --- number of elements in the array wvno
!                wvno --- the first elememnt is wk, the rest elements are
!                        the real roots of the equation um*tahn(um*h)= -v
!                wk --- positive root of the dispersion equation k*tahn(k*h)= v
!       output:  dsa0 --- derivative value of first expansion coefficient by series
!                        representation (m=0)
!
!  contributors list:
!         yingyi liu
!         to be continued...
!

    pure real(rk) function dsa0(r, pa, wk)
        real(rk), intent(in) :: r, pa, wk
        integer n, p
        real(rk) psi, rm, en, sgn, facn

        n = 1
        rm = 100._rk
        dsa0 = 0.0_rk
        sgn = -1.0_rk
        facn = 1.0_rk

! loop: the summation stops at the condition when rm meets the tolerance
!
        do while (dabs(rm) > eps)

            call enxa(n, (r/pa)**2, en)

            psi = 0.0_rk
            do p = 1, n
                psi = psi + 1.0_rk/p
            end do

! calculate the value of each term in the infinite series expansion
!
            rm = 2.0_rk*sgn*(r/pa)**(2.0_rk*n - 1.0_rk)/(facn*pa) &
                 + (wk*pa/2.0_rk)**(2.0_rk*n)*en*2.0_rk*r/(pa**2*facn) - 2.0_rk*sgn*n*wk &
                 *(wk*r/2.0_rk)**(2.0_rk*n - 1.0_rk)*(dlog(wk*r/2.0_rk) + ga - psi)/facn**2 &
                 - 2.0_rk*sgn*(wk*r/2.0_rk)**(2.0_rk*n)/(r*facn**2)

            dsa0 = dsa0 + rm

            n = n + 1
            sgn = -sgn   ! sign function depends on the power n
            facn = facn*n   ! factorial of the numbers from 1 to n

        end do

        dsa0 = dsa0/pi4

    end function dsa0

!  purpose: this program computes the series representations of expansion
!                coefficients (m>=1,m=1,2,3,4...)
!
!  license:
!
!    this routine is part of FINGREEN3D.
!
!    FINGREEN3D is a free software package: you can redistribute it and/or modify it
!    under the terms of the gnu lesser general public license as published by the
!    free software foundation, either version 3 of the license, or (at your option)
!    any later version.
!
!    you should have received a copy of the gnu general public license (either v3
!    or later), along with FINGREEN3D. if not, see <http://www.gnu.org/licenses/>.
!
!  code original author:
!
!    yingyi liu on sept 7, 2013
!
!  modified on:
!
!    july 17, 2018
!
!  reference:
!
!    yingyi liu et al. a reliable open-source package for performance evaluation of
!    floating renewable energy systems in coastal and offshore regions. energy conversion
!    and management, 2018
!
!  parameters:
!       input:   r --- horizontal distance of the source point and the field point
!                pa --- a key parameter determing convergence, i.e., pa= a*h
!                a ---  a key parameter related with convergence
!                m --- number of elements in the array wvno
!                wvno --- the first elememnt is wk, the rest elements are
!                        the real roots of the equation um*tahn(um*h)= -v
!                wk --- positive root of the dispersion equation k*tahn(k*h)= v
!       output:  sam --- value of first expansion coefficient by series representation
!                       (m>=1,m=1,2,3,4...)
!
!  contributors list:
!         yingyi liu
!         to be continued...
!

    pure real(rk) function sam(r, pa, wvn)
        real(rk), intent(in) :: r, pa, wvn
        integer n
        real(rk) rm, en, sgn, facn

        n = 0
        rm = 100._rk
        sam = 0.0_rk
        sgn = 1.0_rk
        facn = 1.0_rk

! loop: the summation stops at the condition when rm meets the tolerance
!
        do while (dabs(rm) > eps)

            call enxa(n + 1, (wvn*pa/2.0_rk)**2, en)

            rm = sgn*(r/pa)**(2.0_rk*n)*en/facn

            sam = sam + rm

            n = n + 1
            sgn = -sgn   ! sign function depends on the power n
            facn = facn*n   ! factorial of the numbers from 1 to n

        end do

        sam = sam/pi4

    end function sam

!  purpose: this program computes the series representations of expansion
!                coefficients with repect to r in linton's method (m>=1,m=1,2,3,4...)
!
!  license:
!
!    this routine is part of FINGREEN3D.
!
!    FINGREEN3D is a free software package: you can redistribute it and/or modify it
!    under the terms of the gnu lesser general public license as published by the
!    free software foundation, either version 3 of the license, or (at your option)
!    any later version.
!
!    you should have received a copy of the gnu general public license (either v3
!    or later), along with FINGREEN3D. if not, see <http://www.gnu.org/licenses/>.
!
!  code original author:
!
!    yingyi liu on sept 7, 2013
!
!  modified on:
!
!    july 17, 2018
!
!  reference:
!
!    yingyi liu et al. a reliable open-source package for performance evaluation of
!    floating renewable energy systems in coastal and offshore regions. energy conversion
!    and management, 2018
!
!  parameters:
!       input:   r --- horizontal distance of the source point and the field point
!                pa --- a key parameter determing convergence, i.e., pa= a*h
!                a ---  a key parameter related with convergence
!                m --- number of elements in the array wvno
!                wvno --- the first elememnt is wk, the rest elements are
!                        the real roots of the equation um*tahn(um*h)= -v
!                wk --- positive root of the dispersion equation k*tahn(k*h)= v
!       output:  dsam --- derivative value of first expansion coefficient by series
!                        representation (m>=1,m=1,2,3,4...)
!
!  contributors list:
!         yingyi liu
!         to be continued...
!

    pure real(rk) function dsam(r, pa, wvn)
        real(rk), intent(in) :: r, pa, wvn
        integer n
        real(rk) rm, en, sgn, facn

        n = 1
        rm = 100._rk
        dsam = 0.0_rk
        sgn = -1.0_rk
        facn = 1.0_rk

! loop: the summation stops at the condition when rm meets the tolerance
!
        do while (dabs(rm) > eps)

            call enxa(n + 1, (wvn*pa/2.0_rk)**2, en)

            rm = 2.0_rk*n*sgn*(r/pa)**(2.0_rk*n - 1.0_rk)*en/(facn*pa)

            dsam = dsam + rm

            n = n + 1
            sgn = -sgn   ! sign function depends on the power n
            facn = facn*n   ! factorial of the numbers from 1 to n

        end do

        dsam = dsam/pi4

    end function dsam

!  purpose: this program computes the expansion coefficients and its
!                derivative in linton's method by integral form.
!
!  license:
!
!    this routine is part of FINGREEN3D.
!
!    FINGREEN3D is a free software package: you can redistribute it and/or modify it
!    under the terms of the gnu lesser general public license as published by the
!    free software foundation, either version 3 of the license, or (at your option)
!    any later version.
!
!    you should have received a copy of the gnu general public license (either v3
!    or later), along with FINGREEN3D. if not, see <http://www.gnu.org/licenses/>.
!
!  code original author:
!
!    yingyi liu on sept 7, 2013
!
!  modified on:
!
!    july 21, 2016
!
!  reference:
!
!    yingyi liu et al. a reliable open-source package for performance evaluation of
!    floating renewable energy systems in coastal and offshore regions. energy conversion
!    and management, 2018
!
!  parameters:
!       input:   r --- horizontal distance of the source point and the field point
!                wvn --- the input wave number
!                no --- the term number to be calculated in linton's expansion
!                t ---  integration variable
!                funtag --- 1: calculate function value; 2: calculate derivative value
!       output:  fun --- value of expansion coefficient or its derivative by integral
!                        representation
!
!  contributors list:
!         yingyi liu
!         to be continued...
!

    pure real(rk) function fun(r, wvn, no, t, funtag)

        integer, intent(in) :: no, funtag
        real(rk), intent(in) :: r, wvn, t

        select case (funtag)

        case (1)
            fun = f1(r, wvn, no, t)

        case (2)
            fun = f2(r, wvn, no, t)

        end select

    end function fun

!  purpose: this program computes the expansion coefficients and its
!                derivative in linton's method by integral form.
!
!  license:
!
!    this routine is part of FINGREEN3D.
!
!    FINGREEN3D is a free software package: you can redistribute it and/or modify it
!    under the terms of the gnu lesser general public license as published by the
!    free software foundation, either version 3 of the license, or (at your option)
!    any later version.
!
!    you should have received a copy of the gnu general public license (either v3
!    or later), along with FINGREEN3D. if not, see <http://www.gnu.org/licenses/>.
!
!  code original author:
!
!    yingyi liu on sept 7, 2013
!
!  modified on:
!
!    july 21, 2016
!
!  reference:
!
!    yingyi liu et al. a reliable open-source package for performance evaluation of
!    floating renewable energy systems in coastal and offshore regions. energy conversion
!    and management, 2018
!
!  parameters:
!       input:   r --- horizontal distance of the source point and the field point
!                wvn --- the input wave number
!                no --- the term number to be calculated in linton's expansion
!                t ---  integration variable
!                funtag --- 1: calculate function value; 2: calculate derivative value
!       output:  fun --- value of expansion coefficient or its derivative by integral
!                        representation
!
!  contributors list:
!         yingyi liu
!         to be continued...
!

    pure real(rk) function funf(v, r, zf, zp, h, t, funtag)

        integer, intent(in) :: funtag
        real(rk), intent(in) :: v, r, zf, zp, h, t

        select case (funtag)

!        integrand of l1
        case (1)
            funf = g1(v, r, zf, zp, h, t)

!        integrand of l1r
        case (2)
            funf = g2(v, r, zf, zp, h, t)

!        integrand of l1z
        case (3)
            funf = g3(v, r, zf, zp, h, t)

        end select

    end function funf

!  purpose: this program computes the value of expansion coefficients
!                in linton's method by integral form.
!
!  license:
!
!    this routine is part of FINGREEN3D.
!
!    FINGREEN3D is a free software package: you can redistribute it and/or modify it
!    under the terms of the gnu lesser general public license as published by the
!    free software foundation, either version 3 of the license, or (at your option)
!    any later version.
!
!    you should have received a copy of the gnu general public license (either v3
!    or later), along with FINGREEN3D. if not, see <http://www.gnu.org/licenses/>.
!
!  code original author:
!
!    yingyi liu on sept 7, 2013
!
!  modified on:
!
!    july 21, 2016
!
!  reference:
!
!    yingyi liu et al. a reliable open-source package for performance evaluation of
!    floating renewable energy systems in coastal and offshore regions. energy conversion
!    and management, 2018
!
!  parameters:
!       input:   r --- horizontal distance of the source point and the field point
!                wvn --- the input wave number
!                no --- the term number to be calculated in linton's expansion
!                t ---  integration variable
!       output:  f1 --- value of expansion coefficient by integral representation
!
!  contributors list:
!         yingyi liu
!         to be continued...
!

    pure real(rk) function f1(r, wvn, no, t)
        integer, intent(in) :: no
        real(rk), intent(in) :: r, wvn, t

        if (dabs(t) < 1.d-8) then

            f1 = 0.0_rk

        else

            if (no == 1) then

                f1 = dexp(-r**2/(4._rk*t) + wvn**2*t)/pi4/t

            else

                f1 = dexp(-r**2/(4._rk*t) - wvn**2*t)/pi4/t

            end if

        end if

    end function f1

!  purpose: this program computes derivative of the expansion coefficients
!                in linton's method by integral form.
!
!  license:
!
!    this routine is part of FINGREEN3D.
!
!    FINGREEN3D is a free software package: you can redistribute it and/or modify it
!    under the terms of the gnu lesser general public license as published by the
!    free software foundation, either version 3 of the license, or (at your option)
!    any later version.
!
!    you should have received a copy of the gnu general public license (either v3
!    or later), along with FINGREEN3D. if not, see <http://www.gnu.org/licenses/>.
!
!  code original author:
!
!    yingyi liu on sept 7, 2013
!
!  modified on:
!
!    july 21, 2016
!
!  reference:
!
!    yingyi liu et al. a reliable open-source package for performance evaluation of
!    floating renewable energy systems in coastal and offshore regions. energy conversion
!    and management, 2018
!
!  parameters:
!       input:   r --- horizontal distance of the source point and the field point
!                wvn --- the input wave number
!                no --- the term number to be calculated in linton's expansion
!                t ---  integration variable
!       output:  f2 --- derivative value of the expansion coefficient by integral
!                        representation
!
!  contributors list:
!         yingyi liu
!         to be continued...
!

    pure real(rk) function f2(r, wvn, no, t)
        integer, intent(in) :: no
        real(rk), intent(in) :: r, wvn, t

        if (dabs(t) < 1.d-8) then

            f2 = 0.0_rk

        else

            if (no == 1) then

                f2 = -r*dexp(-r**2/(4._rk*t) + wvn**2*t)/(8._rk*pi*t**2)

            else

                f2 = -r*dexp(-r**2/(4._rk*t) - wvn**2*t)/(8._rk*pi*t**2)

            end if

        end if

    end function f2

!  purpose: this program computes the value of the survival integral
!                in l1 function in linton's method.
!
!  license:
!
!    this routine is part of FINGREEN3D.
!
!    FINGREEN3D is a free software package: you can redistribute it and/or modify it
!    under the terms of the gnu lesser general public license as published by the
!    free software foundation, either version 3 of the license, or (at your option)
!    any later version.
!
!    you should have received a copy of the gnu general public license (either v3
!    or later), along with FINGREEN3D. if not, see <http://www.gnu.org/licenses/>.
!
!  code original author:
!
!    yingyi liu on sept 7, 2013
!
!  modified on:
!
!    july 21, 2016
!
!  reference:
!
!    yingyi liu et al. a reliable open-source package for performance evaluation of
!    floating renewable energy systems in coastal and offshore regions. energy conversion
!    and management, 2018
!
!  parameters:
!       input:   v --- wave number in deep water , i.e., v= w^2/g
!                r --- horizontal distance of the source point and the field point
!                zf,zp --- z coordinates of the source point and the field point
!                h --- water depth (h>0)
!                t ---  integration variable
!       output:  g1 --- value of the survival integral in l1 in linton's
!                       representation
!
!  contributors list:
!         yingyi liu
!         to be continued...
!

    pure real(rk) function g1(v, r, zf, zp, h, t)
        real(rk), intent(in) :: v, r, zf, zp, h, t
        integer i
        real(rk) sum, val(4)

        if (dabs(t) < 1.d-8) then

            g1 = 0.0_rk

        else

            call chi(zf, zp, h, val)

            sum = 0.0_rk
            do i = 1, 4
                sum = sum + dexp(-v*val(i))*erfcc(val(i)/(2.0_rk*t) - v*t)
            end do

            g1 = dexp(v**2*t**2 - r**2/(4._rk*t**2))*sum/t

        end if

    end function g1

!  purpose: this program computes derivative of the survival integral
!                in l1 function with respect to r in linton's method.
!
!  license:
!
!    this routine is part of FINGREEN3D.
!
!    FINGREEN3D is a free software package: you can redistribute it and/or modify it
!    under the terms of the gnu lesser general public license as published by the
!    free software foundation, either version 3 of the license, or (at your option)
!    any later version.
!
!    you should have received a copy of the gnu general public license (either v3
!    or later), along with FINGREEN3D. if not, see <http://www.gnu.org/licenses/>.
!
!  code original author:
!
!    yingyi liu on sept 7, 2013
!
!  modified on:
!
!    july 21, 2016
!
!  reference:
!
!    yingyi liu et al. a reliable open-source package for performance evaluation of
!    floating renewable energy systems in coastal and offshore regions. energy conversion
!    and management, 2018
!
!  parameters:
!       input:   v --- wave number in deep water , i.e., v= w^2/g
!                r --- horizontal distance of the source point and the field point
!                zf,zp --- z coordinates of the source point and the field point
!                h --- water depth (h>0)
!                t ---  integration variable
!       output:  g2 --- derivative value of the survival integral in l1 with respect
!                       to r in linton's representation
!
!  contributors list:
!         yingyi liu
!         to be continued...
!

    pure real(rk) function g2(v, r, zf, zp, h, t)
        real(rk), intent(in) :: v, r, zf, zp, h, t
        integer i
        real(rk) sum, val(4)

        if (dabs(t) < 1.d-8) then

            g2 = 0.0_rk

        else

            call chi(zf, zp, h, val)

            sum = 0.0_rk
            do i = 1, 4
                sum = sum + dexp(-v*val(i))*erfcc(val(i)/(2.0_rk*t) - v*t)
            end do

            g2 = r*dexp(v**2*t**2 - r**2/(4._rk*t**2))*sum/(2.0_rk*t**3)

        end if

    end function g2

!  purpose: this program computes derivative of the survival integral
!                in l1 function with respect to z in linton's method.
!
!  license:
!
!    this routine is part of FINGREEN3D.
!
!    FINGREEN3D is a free software package: you can redistribute it and/or modify it
!    under the terms of the gnu lesser general public license as published by the
!    free software foundation, either version 3 of the license, or (at your option)
!    any later version.
!
!    you should have received a copy of the gnu general public license (either v3
!    or later), along with FINGREEN3D. if not, see <http://www.gnu.org/licenses/>.
!
!  code original author:
!
!    yingyi liu on sept 7, 2013
!
!  modified on:
!
!    july 21, 2016
!
!  reference:
!
!    yingyi liu et al. a reliable open-source package for performance evaluation of
!    floating renewable energy systems in coastal and offshore regions. energy conversion
!    and management, 2018
!
!  parameters:
!       input:   v --- wave number in deep water , i.e., v= w^2/g
!                r --- horizontal distance of the source point and the field point
!                zf,zp --- z coordinates of the source point and the field point
!                h --- water depth (h>0)
!                t ---  integration variable
!       output:  g3 --- derivative value of the survival integral in l1 with respect
!                       to z in linton's representation
!
!  contributors list:
!         yingyi liu
!         to be continued...
!

    pure real(rk) function g3(v, r, zf, zp, h, t)
        real(rk), intent(in) :: v, r, zf, zp, h, t
        integer i
        real(rk) sum, val(4), dval(4)

        if (dabs(t) < 1.d-8) then

            g3 = 0.0_rk

        else

            call chi(zf, zp, h, val)
            call dchi(dval)

            sum = 0.0_rk
            do i = 1, 4
                sum = sum - v*dval(i)*dexp(-v*val(i))*erfcc(val(i)/(2.0_rk*t) - v*t) &
                      + dexp(-v*val(i))*derfcc(val(i)/(2.0_rk*t) - v*t)*dval(i)/(2.0_rk*t)
            end do

            g3 = dexp(v**2*t**2 - r**2/(4._rk*t**2))*sum/t

        end if

    end function g3
!
!    ===========================================================
!      functions 1/dsqrt(r^2+chi(1,i)^2) and their derivatives
!         code original author: yingyi liu   2013.09.07
!    ===========================================================

    pure subroutine rchi(r, zf, zp, h, ro, dror, droz)
        real(rk), intent(in) :: r, zf, zp, h
        real(rk), intent(out) :: ro(4), dror(4), droz(4)
        integer i
        real(rk) val(4)

        call chi(zf, zp, h, val)

        do i = 1, 4

            ro(i) = dsqrt(r**2 + val(i)**2)
            dror(i) = r/dsqrt(r**2 + val(i)**2)
            droz(i) = val(i)*(-1.0_rk)**i/dsqrt(r**2 + val(i)**2)

        end do

    end subroutine rchi
!
!    ===========================================================
!                        function chi(1,i)
!         code original author: yingyi liu   2013.09.07
!    ===========================================================

    pure subroutine chi(zf, zp, h, val)
        real(rk), intent(in) :: zf, zp, h
        real(rk), intent(out) :: val(4)

        val(1) = -zp - zf
        val(2) = 2.0_rk*h - zp + zf
        val(3) = 2.0_rk*h + zp - zf
        val(4) = 4._rk*h + zp + zf

    end subroutine chi
!
!    ===========================================================
!       the derivative of function chi(1,i) wih respect to z
!         code original author: yingyi liu   2013.09.07
!    ===========================================================

    pure subroutine dchi(dval)
        real(rk), intent(out) :: dval(4)

        dval(1) = -1.0_rk
        dval(2) = 1.0_rk
        dval(3) = -1.0_rk
        dval(4) = 1.0_rk

    end subroutine dchi

!
!    ===========================================================
!                         function nm
!         code original author: yingyi liu   2013.09.07
!    ===========================================================

    pure function nm(h, wvn)
        real(rk), intent(in) :: h
        complex(rk), intent(in) :: wvn
        complex(rk) :: nm

        nm = h*(1.0_rk + cdsin(2.0_rk*wvn*h)/(4._rk*wvn*h))

    end function nm
!
!    ===========================================================
!          derivative of the complementary error function
!         code original author: yingyi liu   2013.09.07
!    ===========================================================

    pure function derfcc(x)
        real(rk), intent(in) :: x
        real(rk) derfcc

        derfcc = -2.0_rk*dexp(-x**2.0_rk)/sqpi

    end function derfcc

!
!    ===========================================================
!             the complementary error function erfc(x)
!         code original author: yingyi liu   2013.09.07
!    ===========================================================

    pure function erfcc(x)
        real(rk), intent(in) :: x
        real(rk) erfcc, err

        call error(x, err)
        erfcc = 1.0_rk - err

    end function erfcc

!  purpose: this program computes the 1d integral value of an external
!                integrand function in a finite interval (a,b) using standard
!                kronrod extension of gauss-lengendre rule,  here using the
!                commonly used standard pair (g7,k15).
!
!  license:
!
!    this routine is part of FINGREEN3D.
!
!    FINGREEN3D is a free software package: you can redistribute it and/or modify it
!    under the terms of the gnu lesser general public license as published by the
!    free software foundation, either version 3 of the license, or (at your option)
!    any later version.
!
!    you should have received a copy of the gnu general public license (either v3
!    or later), along with FINGREEN3D. if not, see <http://www.gnu.org/licenses/>.
!
!  code original author:
!
!    yingyi liu on june 2, 2012
!
!  modified on:
!
!    july 21, 2016
!
!  reference:
!
!    yingyi liu et al. a reliable open-source package for performance evaluation of
!    floating renewable energy systems in coastal and offshore regions. energy conversion
!    and management, 2018
!
!    yingyi liu, ying gou, bin teng, shigeo yoshida
!    an extremely efficient boundary element method for wave interaction with long
!    cylindrical structures based on free-surface green's functions.
!    computation, 4(2016), 36
!
!  parameters:
!       input:   v,r,zf,zp,h --- variables of the integrated function, should be
!                        changed by the user when he/she uses this program separately
!                        from this library for other purpose.
!                a,b --- the two end-points of the interval (a,b).
!                funtag --- a flag for indicating which function components to be integrated.
!  ext. routine: fun --- the external subroutine of the function integrand to be integrated,
!                        this should be specified by the user when he/she uses this program
!                        separately from this library for other purpose.
!       output:  gk_int --- the integral value calculated by gauss-kronrod rule,
!                           which is the final output and more accurate.
!                err --- the error estimate of this integration within interval (a,b).
!
!  passed variables:
!                xgq(7),wgq(7) --- the abscissas and weights for the gauss rule.
!                xkq(15),wkq(15) --- the abscissas and weights for the gauss-kronrod rule.
!                g_int --- the integral value calculated by gauss rule.
!                stf(7) ---  the stored integrand function value calculated by gauss rule,
!                            which can be reused by the gauss-kronrod rule.
!
!  contributors list:
!         yingyi liu
!         to be continued...
!

    pure subroutine gk_integf(v, r, zf, zp, h, a, b, gk_int, err, funtag)
        integer, intent(in) :: funtag
        real(rk), intent(in) :: v, r, zf, zp, h, a, b
        real(rk), intent(out) :: gk_int, err
        integer i
        real(rk) stf(7), g_int

! ** matrix xgq store the gauss-legendre sampling points(7)
!
        real(rk), parameter :: xgq(7) = [-0.949107912342759_rk, -0.741531185599394_rk, &
                                         -0.405845151377397_rk, 0.000000000000000_rk, &
                                         +0.405845151377397_rk, 0.741531185599394_rk, &
                                         +0.949107912342759_rk]
!
! ** matrix wgq store the gauss-legendre weighting factors
!
        real(rk), parameter :: wgq(7) = [0.129484966168870_rk, 0.279705391489277_rk, &
                                         0.381830050505119_rk, 0.417959183673469_rk, &
                                         0.381830050505119_rk, 0.279705391489277_rk, &
                                         0.129484966168870_rk]
!
! ** matrix xkq store the kronrod sampling points(15)
!
        real(rk), parameter :: xkq(15) = [-0.991455371120813_rk, -0.949107912342759_rk, &
                                          -0.864864423359769_rk, -0.741531185599394_rk, &
                                          -0.586087235467691_rk, -0.405845151377397_rk, &
                                          -0.207784955007898_rk, 0.000000000000000_rk, &
                                          +0.207784955007898_rk, 0.405845151377397_rk, &
                                          +0.586087235467691_rk, 0.741531185599394_rk, &
                                          +0.864864423359769_rk, 0.949107912342759_rk, &
                                          +0.991455371120813_rk]
!
! ** matrix wkq store the weighting factors for the kronrod
!    sampling points
!
        real(rk), parameter :: wkq(15) = [0.022935322010529_rk, 0.063092092629979_rk, &
                                          0.104790010322250_rk, 0.140653259715525_rk, &
                                          0.169004726639267_rk, 0.190350578064785_rk, &
                                          0.204432940075298_rk, 0.209482141084728_rk, &
                                          0.204432940075298_rk, 0.190350578064785_rk, &
                                          0.169004726639267_rk, 0.140653259715525_rk, &
                                          0.104790010322250_rk, 0.063092092629979_rk, &
                                          0.022935322010529_rk]

        g_int = 0.0_rk
        do i = 1, 7
            stf(i) = funf(v, r, zf, zp, h, (b + a + (b - a)*xgq(i))/2.0_rk, funtag)
            g_int = g_int + wgq(i)*(b - a)*stf(i)/2.0_rk
        end do

        gk_int = 0.0_rk
        do i = 1, 15
            if (mod(i, 2) == 0) then
                gk_int = gk_int + wkq(i)*(b - a)*stf(i/2)/2.0_rk
            else
                gk_int = gk_int + &
                         wkq(i)*(b - a)*funf(v, r, zf, zp, h, 1/2.0_rk*(b + a + (b - a)*xkq(i)), funtag)/2.0_rk
            end if
        end do
        err = (200.0_rk*dabs(gk_int - g_int))**(1.5_rk)

    end subroutine gk_integf

!  purpose: this program using an adaptive quadrature method to compute the
!                integral value of an external integrand function in a finite
!                interval (a,b) by automatically deviding the interval into halves
!                and continuously calling the gauss-kronrod subroutine in order
!                to finally meet the requested accuracy eps.
!
!  license:
!
!    this routine is part of FINGREEN3D.
!
!    FINGREEN3D is a free software package: you can redistribute it and/or modify it
!    under the terms of the gnu lesser general public license as published by the
!    free software foundation, either version 3 of the license, or (at your option)
!    any later version.
!
!    you should have received a copy of the gnu general public license (either v3
!    or later), along with FINGREEN3D. if not, see <http://www.gnu.org/licenses/>.
!
!  code original author:
!
!    yingyi liu on june 2, 2012
!
!  modified on:
!
!    july 21, 2016
!
!  reference:
!
!    yingyi liu et al. a reliable open-source package for performance evaluation of
!    floating renewable energy systems in coastal and offshore regions. energy conversion
!    and management, 2018
!
!    yingyi liu, ying gou, bin teng, shigeo yoshida
!    an extremely efficient boundary element method for wave interaction with long
!    cylindrical structures based on free-surface green's functions.
!    computation, 4(2016), 36
!
!  parameters:
!       input:   v,r,zf,zp,h --- variables of the integrated function, should be
!                        changed by the user when he/she uses this program separately
!                        from this library for other purpose.
!                a,b --- the two end-points of the interval (a,b).
!                eps --- the requested tolerance for the integration.
!                funtag --- a flag for indicating which function components to be integrated.
!  ext. routine: gk_integ --- the external subroutine of gauss-kronrod rule, this should be
!                        modified a little by the user when he/she uses this program
!                        separately from this library for other purpose.
!       output:  ans --- the result of integration.
!
!  passed variables:
!                gk_int --- the integral value calculated by subroutine gk_integ.
!                err ---  the error estimate calculated by subroutine gk_integ.
!
!  contributors list:
!         yingyi liu
!         to be continued...
!

    pure recursive real(rk) function aquadf(v, r, zf, zp, h, a, b, eps, funtag) result(ans)
        real(rk), intent(in) :: v, r, zf, zp, h, a, b, eps
        integer, intent(in) :: funtag
        real(rk) gk_int, err

        ans = 0.0_rk
        call gk_integf(v, r, zf, zp, h, a, b, gk_int, err, funtag)

        if (err > eps .and. dabs(a - b) > eps) then
            ans = ans + aquadf(v, r, zf, zp, h, a, (a + b)/2.0_rk, eps, funtag) &
                  + aquadf(v, r, zf, zp, h, (a + b)/2.0_rk, b, eps, funtag)
        else
            ans = ans + gk_int
        end if

    end function aquadf

!  purpose: this program computes the 1d integral value of an external
!                integrand function in a finite interval (a,b) using standard
!                kronrod extension of gauss-lengendre rule,  here using the
!                commonly used standard pair (g7,k15).
!
!  license:
!
!    this routine is part of FINGREEN3D.
!
!    FINGREEN3D is a free software package: you can redistribute it and/or modify it
!    under the terms of the gnu lesser general public license as published by the
!    free software foundation, either version 3 of the license, or (at your option)
!    any later version.
!
!    you should have received a copy of the gnu general public license (either v3
!    or later), along with FINGREEN3D. if not, see <http://www.gnu.org/licenses/>.
!
!  code original author:
!
!    yingyi liu on june 2, 2012
!
!  modified on:
!
!    july 21, 2016
!
!  reference:
!
!    yingyi liu et al. a reliable open-source package for performance evaluation of
!    floating renewable energy systems in coastal and offshore regions. energy conversion
!    and management, 2018
!
!    yingyi liu, ying gou, bin teng, shigeo yoshida
!    an extremely efficient boundary element method for wave interaction with long
!    cylindrical structures based on free-surface green's functions.
!    computation, 4(2016), 36
!
!  parameters:
!       input:   r,wvn,no --- variables of the integrated function, should be
!                        changed by the user when he/she uses this program separately
!                        from this library for other purpose.
!                a,b --- the two end-points of the interval (a,b).
!                funtag --- a flag for indicating which function components to be integrated.
!  ext. routine: fun --- the external subroutine of the function integrand to be integrated,
!                        this should be specified by the user when he/she uses this program
!                        separately from this library for other purpose.
!       output:  gk_int --- the integral value calculated by gauss-kronrod rule,
!                           which is the final output and more accurate.
!                err --- the error estimate of this integration within interval (a,b).
!
!  passed variables:
!                xgq(7),wgq(7) --- the abscissas and weights for the gauss rule.
!                xkq(15),wkq(15) --- the abscissas and weights for the gauss-kronrod rule.
!                g_int --- the integral value calculated by gauss rule.
!                stf(7) ---  the stored integrand function value calculated by gauss rule,
!                            which can be reused by the gauss-kronrod rule.
!
!  contributors list:
!         yingyi liu
!         to be continued...
!

    pure subroutine gk_integ(r, wvn, no, a, b, gk_int, err, funtag)
        integer, intent(in) :: no, funtag
        real(rk), intent(in) :: r, wvn, a, b
        real(rk), intent(out) :: gk_int, err
        integer i
        real(rk) stf(7), g_int

! ** matrix xgq store the gauss-legendre sampling points(7)
!
        real(rk), parameter :: xgq(7) = [-0.949107912342759_rk, -0.741531185599394_rk, &
                                         -0.405845151377397_rk, 0.000000000000000_rk, &
                                         +0.405845151377397_rk, 0.741531185599394_rk, &
                                         +0.949107912342759_rk]
!
! ** matrix wgq store the gauss-legendre weighting factors
!
        real(rk), parameter :: wgq(7) = [0.129484966168870_rk, 0.279705391489277_rk, &
                                         0.381830050505119_rk, 0.417959183673469_rk, &
                                         0.381830050505119_rk, 0.279705391489277_rk, &
                                         0.129484966168870_rk]
!
! ** matrix xkq store the kronrod sampling points(15)
!
        real(rk), parameter :: xkq(15) = [-0.991455371120813_rk, -0.949107912342759_rk, &
                                          -0.864864423359769_rk, -0.741531185599394_rk, &
                                          -0.586087235467691_rk, -0.405845151377397_rk, &
                                          -0.207784955007898_rk, 0.000000000000000_rk, &
                                          +0.207784955007898_rk, 0.405845151377397_rk, &
                                          +0.586087235467691_rk, 0.741531185599394_rk, &
                                          +0.864864423359769_rk, 0.949107912342759_rk, &
                                          +0.991455371120813_rk]
!
! ** matrix wkq store the weighting factors for the kronrod
!    sampling points
!
        real(rk), parameter :: wkq(15) = [0.022935322010529_rk, 0.063092092629979_rk, &
                                          0.104790010322250_rk, 0.140653259715525_rk, &
                                          0.169004726639267_rk, 0.190350578064785_rk, &
                                          0.204432940075298_rk, 0.209482141084728_rk, &
                                          0.204432940075298_rk, 0.190350578064785_rk, &
                                          0.169004726639267_rk, 0.140653259715525_rk, &
                                          0.104790010322250_rk, 0.063092092629979_rk, &
                                          0.022935322010529_rk]

        g_int = 0.0_rk
        do i = 1, 7
            stf(i) = fun(r, wvn, no, (b + a + (b - a)*xgq(i))/2.0_rk, funtag)
            g_int = g_int + wgq(i)*(b - a)*stf(i)/2.0_rk
        end do

        gk_int = 0.0_rk
        do i = 1, 15
            if (mod(i, 2) == 0) then
                gk_int = gk_int + wkq(i)*(b - a)*stf(i/2)/2.0_rk
            else
                gk_int = gk_int + wkq(i)*(b - a)*fun(r, wvn, no, (b + a + (b - a)*xkq(i))/2.0_rk, funtag)/2.0_rk
            end if
        end do

        err = (200.0_rk*dabs(gk_int - g_int))**(1.5_rk)

    end subroutine gk_integ

!  purpose: this program using an adaptive quadrature method to compute the
!                integral value of an external integrand function in a finite
!                interval (a,b) by automatically deviding the interval into halves
!                and continuously calling the gauss-kronrod subroutine in order
!                to finally meet the requested accuracy eps.
!
!  license:
!
!    this routine is part of FINGREEN3D.
!
!    FINGREEN3D is a free software package: you can redistribute it and/or modify it
!    under the terms of the gnu lesser general public license as published by the
!    free software foundation, either version 3 of the license, or (at your option)
!    any later version.
!
!    you should have received a copy of the gnu general public license (either v3
!    or later), along with FINGREEN3D. if not, see <http://www.gnu.org/licenses/>.
!
!  code original author:
!
!    yingyi liu on june 2, 2012
!
!  modified on:
!
!    july 21, 2016
!
!  reference:
!
!    yingyi liu et al. a reliable open-source package for performance evaluation of
!    floating renewable energy systems in coastal and offshore regions. energy conversion
!    and management, 2018
!
!    yingyi liu, ying gou, bin teng, shigeo yoshida
!    an extremely efficient boundary element method for wave interaction with long
!    cylindrical structures based on free-surface green's functions.
!    computation, 4(2016), 36
!
!  parameters:
!       input:   r,wvn,no --- variables of the integrated function, should be
!                        changed by the user when he/she uses this program separately
!                        from this library for other purpose.
!                a,b --- the two end-points of the interval (a,b).
!                eps --- the requested tolerance for the integration.
!                funtag --- a flag for indicating which function components to be integrated.
!  ext. routine: gk_integ --- the external subroutine of gauss-kronrod rule, this should be
!                        modified a little by the user when he/she uses this program
!                        separately from this library for other purpose.
!       output:  ans --- the result of integration.
!
!  passed variables:
!                gk_int --- the integral value calculated by subroutine gk_integ.
!                err ---  the error estimate calculated by subroutine gk_integ.
!
!  contributors list:
!         yingyi liu
!         to be continued...
!

    pure recursive real(rk) function aquad(r, wvn, no, a, b, eps, funtag) result(ans)
        real(rk), intent(in) :: r, wvn, a, b, eps
        integer, intent(in) :: no, funtag
        real(rk) gk_int, err

        ans = 0.0_rk
        call gk_integ(r, wvn, no, a, b, gk_int, err, funtag)

        if (err > eps .and. dabs(a - b) > eps) then
            ans = ans + aquad(r, wvn, no, a, (a + b)/2.0_rk, eps, funtag) &
                  + aquad(r, wvn, no, (a + b)/2.0_rk, b, eps, funtag)
        else
            ans = ans + gk_int
        end if

    end function aquad

end module fin_green_3d_m
